Transient regime in non-linear transport through many-level quantum dots 
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We investigate the nonstationary electronic transport in noninteracting nanostructures driven 
by a finite bias and time-dependent signals applied at their contacts to the leads. The systems 
are modelled by a tight-binding Hamiltonian and the transient currents are computed from the 
non-equilibrium Green-Keldysh formalism. The numerical implementation is not restricted to weak 
coupling to the leads and does not imply the wide-band limit assumption for the spectral width of 
the leads. As an application of the method we study in detail the transient behavior and the charge 
dynamics in single and double quantum dots connected to leads by a step-like potential, but the 
method allows as well the consideration of non-periodic potentials or short pulses. We show that 
when the higher energy levels of the isolated system are located within the bias window of the leads 
the transient current approaches the steady state in a non-oscillatory smooth fashion. At moderate 
coupling to the leads and fixed bias the transient acquires a step-like structure, the length of the 
steps increasing with the system size. The number of levels inside a finite bias window can be tuned 
by a constant gate potential. We find also that the transient behavior depends on the specific way 
of coupling the leads to the mesoscopic system. 

PACS numbers: 73.23.Hk, 85.35.Ds, 85.35.Be, 73.21.La 
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I. INTRODUCTION 

The dynamics of conduction electrons in open nanos- 
tructures modulated by time-dependent signals is an out- 
standing problem in quantum transport theory. Exten- 
sive experimental and theoretical work has been done es- 
pecially on quantum pumpingi*2iMiSi& (a detailed bibli- 
ography can be found in the recent review^) and photon- 
assisted tunneling (see Ref. [8| and references therein). 
In these phenomena one is interested in measuring or 
computing the current response of a mesoscopic system 
driven by time-dependent periodic signals applied either 
on the system or on the attached leads. In particular, an 
unbiased system subjected to two periodic potentials dif- 
fering by a phase lag generates a nonvanishing pumped 
current, provided one averages over the relevant period. 4 
If the signal frequency is small the pumping is adiabatic 
and can be described by a 'frozen' scattering matrix as 
is rigorously shown in Ref. ||| A photon-assisted tun- 
neling implies instead high frequencies and the measured 
current displays satellite peaks due to the additional side- 
bands^ 

From the theoretical point of view, any calculation of 
the current starts from defining the initial equilibrium 
state of the system and the perturbation that drives it. 
One way is to start with the connected system in the ab- 
sence of the bias, and then to apply the bias adiabatically 
performing linear response calculations for the steady 
state current. An alternative picture was proposed by 
Caroli— that takes as the equilibrium state the state of 
the decoupled system with the bias already imposed on 
the leads. The perturbation here is instead the coupling 
to the leads that is usually adiabatically switched in the 
remote past and eventually reaches its full magnitude at 
t = Q. If one assumes that a steady state is achieved, 
the Green functions depend only on time differences and 



the Keldsyh formalism gives the corresponding current in 
terms of Fourier transformed quantities. Both pictures 
were shown to be useful in capturing and explaining im- 
portant effects. 

As for an ac signal, it can be included in the Green- 
Keldysh approach as a time-dependent global shift of the 
spectrum of the leadsJ^ Note that the occupation prob- 
ability (i.e. the Fermi function) of the leads stays time- 
independent. Therefore this procedure assumes somehow 
that an ac signal is applied as well adiabatically. In the 
quantum pumping calculations the pumping potential is 
applied to the system in a steady-state, be it subjected to 
a finite bias or not. Finally the current (transient, time- 
averaged or steady-state) is computed from the Keldysh- 
Green function formalism^ 

Here we aim to get some insight into a related topic: 
the calculation of the transient current through a quan- 
tum dot (QD) whose coupling to the leads is time de- 
pendent while the bias applied on the leads is constant. 
There are considerably fewer theoretical results on this 
issue (see the references below) and we were motivated 
also by recent increasing interest in using suitable elec- 
tric pulses to investigate relaxation processes in quantum 
dots by using pump and probe measurements or tran- 
sient current spectroscopy^ As in the well known case 
of a turnstile pumpi^ these techniques imply oscillating 
tunnel barriers so that the transport formalism should 
deal with the nonlocal time-dependent coupling between 
the leads and the system. 

The problem we want to look at is defined as follows: 
i) The system is disconnected at any time t < 0, and the 
leads are submitted to a constant bias which is included 
through the difference between the chemical potentials of 
the leads, ii) At t = the leads are suddenly plugged to 
the system. Physically this means that the tunneling bar- 
riers between the leads and the system are set to be very 
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high at t < and drop suddenly at t — to an interme- 
diate value that allows charge transfer across the system. 
The simplest case of a constant barrier height at any 
t > opens already the problem of the existence of non- 
equilibrium steady-states in the long time limit. In gen- 
eral and at a rigorous level, one can prove the existence 
of such states when t — > oo^ and moreover, a Landauer 
formula was shown to hold for the steady-state current ^ 
The method we developed is able to check the passage 
from transient behavior to steady state regime for specific 
systems, like many-level one- and two-dimensional quan- 
tum dots. As we shall see, the onset of the steady state 
for a given system depends on the its structure (number 
of levels), on the measurement setup (the strength of the 
coupling to the leads and the location of the contacts), 
and also on external parameters like gate potentials. 

Although, in the present work the numerical simula- 
tions are restricted to the step- like coupling to the leads, 
our model allows the consideration of more general time- 
dependent potentials between the leads and the central 
region. In particular we can investigate the response of 
a system to nonlocal time-dependent perturbations that 
can be switched on and off individually. 

Recently there have been several theoretical ap- 
proaches to the transient regime. In Ref. [19| the time- 
dependent density functional was used to compute the 
transient current in a onc-dimcnsional system submitted 
to a finite bias applied on the leads. The coupling term 
does not depend on time and the system is in an equi- 
librium state in the absence of the bias. Starting from 
this state the Kohn-Sham equation is used to calculate 
the response of the system to the external bias. The 
same techniques allow the calculation of time-averaged 
current in one-dimensional quantum pumps ^ On the 
other hand, Maciejko et al^ have computed within the 
Keldysh framework the response of a single-site dot for a 
step-like or periodic signal applied to the leads, beyond 
the wide band limit. In their approach the computation 
of the time-dependent Green functions of the perturbed 
system uses steady-state Green functions of the coupled 
and biased system that have to be provided from density 
functional theory. We believe that theoretical results re- 
garding the transient regime induced by time-dependent 
couplings of many-level systems would complement these 
results. 

The content of the paper is organized as follows: Sec- 
tion II presents the model and the theoretical tools we use 
to compute the transient current. We rely essentially on 
the non-equilibrium Green-Keldysh formalism. However, 
in contrast to most of the previous studies we allow a 
complex structure for the central region coupled to leads 
(i.e. there is more than one single localized level and the 
system can be as complicated as we want: a single dot, a 
double dot or an Aharanov-Bohm interferometer). Also 
we go beyond the wide-band limit approximation and we 
solve exactly the integral Dyson equation for the retarded 
Green function of the coupled central region by a suitable 
numerical procedure. In doing so we take into account 



all the scattering processes between the leads and the 
sample. 

Although the electron-electron interaction could pre- 
sumably play an important role in the transient behavior 
and the formalism we use allows the inclusion of Coulomb 
terms in the Hamiltonian, we do not take it into ac- 
count in this work. As is well known, the problem of the 
Coulomb interaction in the Keldysh approach is mainly 
technical and implies suitable approximation schemes for 
the interaction self-energy. We postpone this issue for fu- 
ture work. Section III gives extensive discussion of the 
numerical simulations for single and double dots. Section 
IV concludes the paper. 



II. FORMALISM 

The systems we study in this work have a typical trans- 
port configuration: a central region (S) coupled to two 
semiinfinite leads (a and (3) via a tunneling term (see 
Fig. 1 for a schematical representation) . We shall use a 
tight-binding (TB) description of the Hamiltonian which 
has the following form: 



H{t) = H s + H L + H T (t), 



(1) 



where H$ describes the system, Hi, the semiinfinite leads 
and Hx{t) is the time-dependent tunneling term: 



H T{t)= J2 V *™{t){c\d m + h.c) 



(2) 



Here Ci and cj denote the annihilation/creation opera- 




v„(t) 




v R (t) 



FIG. 1: Schematic picture of the system. The step-like po- 
tential is applied between the leads and the central region S 
at t = 0. 

tors on the z-th site of the lead 7. Similarly d m and d} n 
is the pair of operators corresponding to the TO-th site 
from the central region S. Vi rn (t) is the time-dependent 
hopping coefficient between the i-th site of the lead 7 
and the m-th site of the central region. We take here 
a nearest-neighbor coupling so the double sums in the 
above expression contain only pairs of sites from the leads 
endpoint and the corresponding contact region from the 
central system: 



V im (t) 



VJt) 



ifi,m nearest neighbors , , 
otherwise ^ ' 
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In this work we consider a steplike potential, i.e. V 1 (t) = 
V 1 if t > and zero otherwise. H$ has a usual tight- 
binding form 

N 

H S =Y1 ( £m + V 9)4n d m + tmncQndn. (4) 

m=l (m,n) 

Here t mn are hopping terms, (m, n) denotes nearest- 
neighbor summation over the system sites. e m is the on- 
site energy and the diagonal term V g simulates a plunger 
gate potential applied on the system. N is the number 
of sites in the central region. The spectral width of the 
tight-binding lead is as usual w :— [— 2t£,2ti], where tz, 
is the hopping energy on leads (we take the same hopping 
constant on every lead) . In the numerical calculation we 
choose tL such that it covers entirely the spectrum of the 
central region but we do not assume it to be infinite, as 
it is done in the wide-band limit approximation. 

The central problem in electronic transport is to com- 
pute the statistical average of the time-dependent current 
operator in a given lead (say a) J a (t) = Tr{p(t)j a (t)} 
using the statistical operator p{t). Notice that the time- 
dependence of the current operator appears only because 
of the time-dependent coupling. Denoting by {i a } the 
endpoint sites of the lead a which are coupled to the 
sites {m a } of the central region and by M the number of 
sites in the transverse direction of the lead, the current 
operator j a has the form (we take the electron charge as 
— e and e > 0): 

M 

**(*) = ^ J2 J2 v ^At)(4j ma -dl a c t J. (5) 

i a =l m a eC a 

Since the statistical operator p(t) of the coupled system 
is not easy to compute it is useful to move the time- 
dependence entirely to the current operator by writing 
p(t) in terms of the equilibrium statistical operator po of 



the disconnected system. In general, the coupling to the 
leads is established at a given instant to so that p(t) = po 
for t < <o- Then using the unitary evolution U of the 
full Hamiltonian in the interaction picture w.r.t. the 
unperturbed one the solution of the quantum Liouville 
equation is given by 

p{t) = e-^ Hs+H ^U(t, t )p(t )U(t, t ye^ Hs+H ^ (6) 

Then it can be shown (see e.g. Ref. [2~2~| ) that: 

Ut) = Tr{poTcie~ lI c ds "^3 a (t))}, (7) 

where Tc is the ordering operator on the Schwinger- 
Keldysh contour C that runs from to to t and back to to. 
We remind the reader that in the case of adiabatic cou- 
pling the statistical operator becomes time-independent 
only in the remote past to — > — oo. Both the coupling 
and current operators are written in the interaction pic- 
ture. Using the definitions of the lesser Green functions 
in terms of the Heisenberg operators: 

G< aia (tJ) = i(ct(t')d ma (t)) 
Gf ama M = i(dljt')c ioi (t)), (8) 

it follows that the current is given by a simpler relation: 
2 M 

ai^E ^ MV iam At)G< aia (t,t)). (9) 

At this point the standard Keldysh formalism requires 
the application of the so called Langreth rulesi^ in order 
to express the lesser Green function G^ oio in terms of the 
Green functions of the central region in the presence of 
the leads and the Green functions of the isolated 

semiinfinite lead gf'f- The latter can be analytically 
computed: 



gf a ja (t, t') = i6{t' - t) V Xp(i a )x P (ja) / " ' " dEp(E - E p ) e - iE ^ (10) 

p=l J-2t L +E p 
M r 2t L +E p 

9tj a W) = i^Xp(i a )Xv(Ja) / dEp{E-E p )e- iE ^f a {E). (11) 

p— 1 J —2t.L-\-E p 

i 



In the above equations E p = Hl cos(ptt/(M + 1)) is the 
energy of the transverse channel p. These channels ap- 
pear due to the width of the leads which in the tight- 
binding description is given by the number of the sites 
M in the transverse direction. More exactly, the many 
channel lead is constructed by taking M semiinfinite ID 
leads and by coupling them through nearest neighbor 



hopping constants. Then Xp(*q) — y m+i s ^ n ^ M+i ) ^ s 
the transversal eigenfunction associated to E p , and p{E) 
is the density of states at the endpoint of a semiinfinite 
one-dimensional lead: 



y/U 2 r - E 2 

p(E)=8(2t L -\E\)^_ . (12) 
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Finally, f a (E) is the Fermi function in the lead a. The 
bias is included in our approach as the difference between 
the two chemical potentials of the leads V = hl — Hr- 

I 



Ja(t) = 



We have introduced the energy and time-dependent 
quantity r^'^ that takes also into account the M chan- 
nels in the lead: 

M 

r%ln a ( E ;t,s)= £ p(E-E p ) Xp (i a ) Xp (j<*)x 

V iama (t)V jana (s). (14) 
A similar formula can be written down for the current Jp 
I 



G R {t,t') 
G<(t,t') 



where G R ' A {t,t') are the retarded and advanced Green 
functions of the isolated central region and Y, R,< are the 
retarded and lesser self-energies. G R (t,t') has a simple 
expression in terms of the discrete spectrum {E\} of the 
central region and its localized eigenfunctions ip\ (clearly 
\=1,..N): 

G R >mn {t,t') = -i6(t - f ) J2 ^H**! 1 -''). 

We emphasize the lower integration limit t = in equa- 
tions (fl6|) and (fl~7|) . This is due to the fact that there is no 
coupling for t < 0. In the adiabatic setup the coupling 
is established in the remote past and one should set a 
lower cutoff in the numerical implementation. However, 
the Dyson equation still contains two coupled integrals. 
The two self-energies above contain the information from 
the leads and are finite rank matrices in the Hubert space 
of the central region S (7 = a, f3): 

7 

££»(*,*') = E K rW5, r ;^(^ t, )^(i')^ m ^™n T (20) 

7 



Plugging all these elements in the current formula one 
gets the main expression that will be numerically imple- 
mented in the next section: 



and therefore one defines as well the net current 

J(t) = J a {t) + J p (t). (15) 

It is important to observe that in contrast to the simple 
case of a single-site system the expressions for the two 
currents imply the Green functions at different contacts. 
The retarded and the lesser Green functions are then to 
be computed from the Dyson and Keldysh equations^: 

(16) 
(17) 

We stress that the indices of the leads' Green function are 
unambiguously determined as the neighbor sites of the 
contact surface C a ■ In the single channel case M = 1 one 
recovers simpler expressions. In particular the retarded 
Green function of the lead can be expressed through the 
Bessel function of the first kind: 

5 W*>*J- 2t L (t-t>) ' [ > 

We point out the difference between the exact form of 
the retarded self-energy and the simple wide-band limit 
expression (which simplifies to S(t — t') up to some con- 
stants). Note that the retarded Green function gives the 
leads' self-energy and is a highly oscillating functions, it 
will turn out in Section III that this behavior has cru- 
cial effects on the transient current. Another difficulty of 
Eq. comes from the quadratic dependence of the self- 
energies on the time-dependent coupling. Clearly this 
prevents any partial Fourier transform trick. 

Given these, our strategy in solving the integral Dyson 
equation relies in transforming it into an algebraic equa- 
tion of the form AX = B where A, X, B are generalized 
complex matrices depending on both spatial and time 
arguments. To this end we first plug the retarded self 
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XL £ 



P— 1 rn Q ,n a GCo 



2t L +E p 



dE 



2tL-\-E v 



dse -iE{s-t) T ^ (E;t,s)(G* (t,s)f a (E) +G< (t,s))). (13) 



G R (t,t')+ / dtiG R (t,ti) / dt 2 X R {t 1 ,t 2 )G R (t 2 ,t') 
Jo Jo 

dttG^tt) f dt 2 Y 1 < {t ll t 2 )G A {t 2 ,t'), 
Jo 

I 
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energy from Eq. (fl9|) in the Dyson equation (fT6|) and dis- 
cretize the time arguments. Note that the variable t 2 is 
defined on a denser grid than the one used for t%. The 
inner time integral is evaluated by a repeated 4-point 
Gauss method, which turned out to be accurate enough 
for the numerical results to be stable when increasing 
the number of integration steps. This procedure allows 
us to write the double integral as a matrix G R A, where 
A is actually a product of G R , T^ R and some diagonal 
matrices containing the Gauss weights needed in the in- 
tegration procedure. Then the adjoint of the general- 
ized retarded Green function G R is simply the solution 

of the algebraic equation (1 — A)*G R = G R . The true 
Green function is recovered by turning back the mixed 
indices of G R . We stress that by solving the equation for 
G R the Dyson equation is solved exactly. Moreover, no 
matrix inversion is required. This is certainly an advan- 
tage in the numerical simulations since it is known that 
matrix inversion is both memory and time-consuming. 
The advanced Green function is computed using the iden- 
tity G^(t,t') — G R (t',t) and the lesser Green function 
is derived from the Keldysh equation. Also, the time- 
dependent occupation number can be computed as: 



N(t) 



Im G 

m£S 



< 

mm 



(22) 



The current in the right lead Jp has a similar expression. 
We note that for a system with many sites one has to deal 
with different contact Green functions besides replacing 
only the Fermi function in the first term in the current 
formula. Moreover, in the transient regime the current 
conservation does not imply the usual identity J a = — J p. 
We shall discuss this feature below. 

In the Keldysh approach to time-dependent transport 
the problem is to extract physical information from the 
two contributions in Eq. (|13p . In the simplest case of a 
single-site and within the wide-band limit it was shown 
that the average current obeys a Landauer-like formula. 
The effects of a step-like or harmonic time-dependent po- 
tentials applied adiabatically on leads were studied both 
in the WBLi£ and beyoncUS. However, to our best knowl- 
edge no transient current calculation for a many-level 
structure beyond the wide-band limit has been performed 
within the Keldysh formalism. 



III. NUMERICAL SIMULATIONS 



lead [—2^,2^]. We take also e = % = 1. The current 
given by Eq. (|13[) can be written as a sum of two contri- 
butions 



Mt) = J^(t) + J<(t), 



(23) 



the '<' and 'R' labelling emphasizing that the corre- 
sponding term contains the lesser and retarded Green 
functions. We shall consider for simplicity only single 
channel leads. 
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FIG. 2: (Color online) (a) The transient current for one site 
coupled to single-channel leads. We present curves for dif- 
ferent values of the coupling strength U and of the bias V. 
(b) The effect of the coupling amplitude U on the occupation 
number. The bias is fixed to V = 1.0, and kT = 0.0001. 



In all the plots the bias, the energy, the hopping con- 
stants on the leads, the coupling strengths and the gate 
potentials will be expressed in terms of the hopping en- 
ergy of the central region to which is chosen as energy 
unit. The current is therefore given in units of etn/Tl 
and the time expressed in units of 1/tjj. Since the spec- 
trum of the two-dimensional discrete Laplacian covers 
the range [—4^,4^] we shall take = 2 in order to 
match it to the spectral width of the one-dimensional 



We start this section by discussing the case of a single 
site dot which allows qualitative discussion on the tran- 
sient regime and on the transition towards the steady- 
state. The site is coupled to single-channel leads (i.e. 
N = M = 1). Both leads are coupled suddenly to the 
system with the same strength U, i.e V a — Vp := U. 
The bias is applied symmetrically on leads, i.e. [i a ,f3 = 
/io±eF/2, /^o being the chemical potential of the unbiased 
leads. We observe that for hq = 0.0 the single eigenvalue 
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of the isolated system Eq = 0.0 is located in the middle 
of the bias window W — [fi — eV/2, [i + eV/2]. As we 
shall see later on, the position of the eigenvalues of the 
system within the bias window has important implica- 
tions on the transient current. The free retarded Green 
function of the single site system is simply G§ = —i but 
the full retarded Green function is still given by the in- 
tegral Dyson equation and an analytic solution is not at 
hand. 

Fig. 2(a) gives the transient current for different val- 
ues of the bias V and of the coupling amplitude U and 
reveals that the parameter that controls the shape and 
the amplitude of the oscillation is the coupling strength 
U. At moderate coupling U = 0.75 the steady state (SS) 
is achieved fast but an oscillatory behavior is observed 
at U = 0.95. The case U — 1.15 is beyond the per- 
turbative regime and the steady-state is not achieved in 
the selected time-range. As the bias increases the cur- 
rent saturates for values of V that exceed the spectrum 
of the leads (i.e. for V > 8), emphasizing the non linear 
transport regime. In turn, the bias neither affects the 
amplitude nor the period of the oscillations. This is due 
to the fact that in our model, as in all approaches based 
on the Keldysh formalism, there is no term in the Hamil- 
tonian to describe the voltage drop across the sample, 
the bias being included only via the Fermi functions of 
the leads. 
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FIG. 3: (Color online) The two contributions to the transient 
current in the left lead. J a if) is always positive while the 
lesser contribution J a it) is negative. The chosen values for 
the coupling strength U are given in the figure. The bias 
V = 1.0, and kT = 0.0001. 

In Fig. 2(b) we plot the occupation number of the res- 
onant site N(t) for the same parameters as in Fig. 1(a). 
The behavior w.r.t. U is similar. Few more things worth 
to be noticed: i) In the steady state regime the occupa- 
tion of the site is 1/2; ii) as U increases while the bias 
stays constant the occupation number reaches its max- 
imum value faster and the value corresponding to the 
steady state decreases; hi) Comparing Fig. 2(a) and 2(b) 
it can be seen that there is no clear relation between the 
principal maximum of the current and the one of the oc- 



cupation number; in fact the electrons are accumulating 
in the system even after the current starts to decrease 
towards the steady state. 

Fig. 3 shows the two contributions to the current j£ 
and J < for V — 1.0 and several coupling constants con- 
sidered in Fig. 2. A physical significance of these two 
currents was proposed in Ref. LL2 for the single site case. 
Although both Green functions at the contacts appear- 
ing in the current formula are 'dressed' by the leads' self- 
energy one could view as the current flowing towards 
the sample and J< (t) as the current from to sample to 
the lead a. One notices that the currents have opposite 
signs. Another observation is that the lesser contribu- 
tion is responsible for the total current oscillations since 
saturates quickly. However at small times grows 
faster than J<, leading thus to the fast increase of the 
transient. 

In order to understand the nature of the oscillations in 
the transient current and their dependence on the cou- 
pling strength U it is useful to rewrite the current formula 
Eq. (I13p in a more useful form (since we consider a single 
site system there is only one contact site and the indices 
of the Green functions can be omitted): 



J a (t) = -2U J Im / ds(G H (t,s)F 1 (s,t)+G < (t,s)F 2 (s,t)), 

JO 

where Fi, F2 are two oscillating integrals: 



21 L 



Fi(s,t) = / dEf a (E)p(E)e 

J-2t L 

F 2 (s,t) = / dEp(E)e- iE ^ 
J -it,, 



iE(s-t) 



(25) 



(26) 



One can easily observe that actually F 2 (s,t) can be ex- 
pressed through Bessel function of the first kind: 



F 2 (s,t) 



9{t- s).h{2t L {t~ s)) 
2t L (t-s) 



(27) 



For fixed t, F 2 is an oscillating function of s whose os- 
cillation amplitudes increase with s. F\ does not have a 
simple analitical expression but it has a similar behavior. 

The oscillatory behavior of the current is clearly de- 
cided by the convolution in Eq. (|2"3|) . Besides the os- 
cillations of i*i ans F 2 one expects as well a complex 
behavior of the Green functions. We recall that the 
Dyson equation counts the infinite back-and-forth tun- 
neling processes involving the leads and that the am- 
plitudes of these events are even powers of U. Now, 
the higher order terms in the Dyson equation contain 
multiple integrals of products of the leads' self-energy 
which is highly oscillating (see (|21jl) Therefore if U << 1 
there will be only few low-order significant contributions 
from the complicated lead-sample scattering. The critical 
value U = 1.00 corresponds to the onset of the nonper- 
turbative regime, and the method we use for solving the 
Dyson equation captures as well this situation, taken into 
account all contributions. 
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FIG. 4: (Color online) The imaginary parts of the retarded (a) 
and lesser (b) two-time Green functions of the coupled single 
site. The bias V = 1.0 and the coupling strength U = 1.2. 
The oscillations seen along the 'diagonal' in (b) corresponding 
to almost equal times are responsible for the oscillations of the 
occupation number. 



We give in Fig. 4 the 3D plots of the imaginary parts 
for the retarded and lesser Green functions at coupling 
strength U — 1.20, which leads to oscillations of the tran- 
sient. These are the relevant quantities in the current 
formula since it turns out that the real part of G R and 
of Fi are vanishingly small (not shown). 

One observes that G R (t,s) — for s > t and, more 
interestingly, that G R (t,s) and G < (t,s) exhibit pro- 
nounced oscillations as time varies and reach a limit value 
as s approaches t. In the case of the retarded Green 
function this limit is constant and equals —i, which is 
simply the value of the unperturbed retarded Green func- 
tion. This feature is easy to understand by looking at the 
Dyson equation and noticing that when s — > t the inte- 
gration range of the inner integral shrinks considerably so 
that at almost equal times the perturbed Green function 
resembles the unperturbed one. This argument is not re- 
stricted to the single site case we are discussing. Also, 
since the spectrum of the discrete Laplacian is symmet- 
ric, the unperturbed retarded Green function will always 
be real and therefore the real part of the full Green func- 
tion will always be vanishingly small, as we shall check 
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FIG. 5: (Color online) The imaginary part of the lesser two- 
time Green function for the moderate coupling U = 0.75 (a) 
and strong coupling U — 1.2 (b). (c) The real part of the 
oscillating function F^. The bias V = 1.0. 



numerically in the many site case. In contrast, the limit 
of ImG < as s — > t is not a constant w.r.t. t but shows 
oscillations that disappear as t increases. In the partic- 
ular single-site case the limit value of G < is clearly the 
occupation number of the site, whose oscillations were 
shown already in Fig. 2(b). 

Figs. 5 (a) and (b) show 3D maps of the imaginary part 
of the lesser Green function and emphasize the role of the 
coupling strength on the transient. At moderate coupling 
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B. Many-site case. 
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FIG. 6: (Color online) (a) The transient current in the left 
lead J a (t) for different sizes of the ID central region. The 
number of sites N is indicated in the figure. In (a) the cou- 
pling to the leads is U = 0.75 and in (b) U = 0.50. By 
decreasing U the shoulders in (a) turn to clear steps in (b). 
The bias is fixed to V = 2.0, and kT = 0.0001. 



U = 0.75 one observe small amplitude oscillations except 
for s ~ i, in clear contrast to the case U = 1.20 where 
pronounced oscillations exist even for large time differ- 
ences. Inspecting the real part of the function F2 given 
in Fig. 5 (c) we see that it does not depend on t when 
s ~ t and that this gives the main contribution to the 
integral It is now clear from Figs. 5 (a) and (c) 

that the corresponding current will be nearly stationary 
once the sample is charged (i.e. for t > 0.75), because 
by increasing t the 'off-diagonal' contributions are very 
small (some cancelations being possible as well). When 
U increases the integral will collect instead nonnegligible 
contributions from the entire range (Q,t) and therefore 
the current will oscillate. These observations lead to the 
following statement: The steady state will be achieved 
at instant t s if for any t > t s there are no contributions 
for long-time differences, i.e. when both contact Green 
functions G R (t,s) and G < (t,s) vanish for s > t s . It is 
easy to observe that this condition implies the well known 
criteria for the steady state G(t, s) = G(t — s, 0). 
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FIG. 7: (Color online) The imaginary part of the left contact 
retarded Green function (a) and of the lesser Green func- 
tion for the 4 site system (b). (c) The imaginary part of 
Gii(ti = 10,t2) at different couplings U. At small coupling 
lmG R resembles the retarded Green function of the isolated 
system (the curve corresponding to U = 0.0). 

We consider now the more interesting case where the 
central region has more than one level. Fig. 6(a) empha- 
sizes the qualitative differences between the transients 
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of ID systems with N — 2,4,6,8 sites. A general fea- 
ture is that as the size increases the transient develops a 
'shoulder' which is not met in the single site case. For 
N = 2 a second smaller slope of decrease is noticed for 
t G [1.6 : 2.5]. For N = 4 the system experiences few very 
different regimes before reaching the steady state. It first 
decreases faster up to t ~ 1.25. For a short time a small 
hill develops around t ~ 2 then the decrease continues 
but clearly at a smaller rate. Finally, for t > 3.25 the cur- 
rent approaches the steady state very slowly. A similar 
behavior is observed for N — 6 and N — 8, the main dif- 
ference being that the 'shoulder' is longer and the inter- 
mediate slope is smaller. The patterns described above 
suggest that there are some intermediate regimes, some 
of them being characterized by a rather stable current. 
Fig. 6(b) emphasizes that at lower coupling U = 0.50 the 
transient is even smoother and for N — 4, 6 and 8 one 
notices the formation of clear steps. 

When the coupling strength U is increased the tran- 
sient shows oscillations but they are fewer than in the 
single site (not shown). Since the two oscillatory inte- 
grals over energy in the current formula and both self- 
energies do not depend on the number of sites in the 
system the above size effects should be explained only by 
the behavior of the contact Green functions. We show in 
Fig. 7(a) and (b) the imaginary parts of the the contact 
Green functions G n ' of the 4 site system for a strong 
coupling U = 0.25 (it turns out again that the real part 
of Gii is vanishingly small). Comparing with Fig. 4 it 
is obvious that for the 4 site dot the Green functions 
have a more regular behavior and in particular the oc- 
cupation number of the contact site 1 shows milder os- 
cillations than the occupation number of the single site. 
Fig. 7(c) gives the ImGf 1 (ii = 10, ti) as a function of t 2 
for different couplings and reveals that at weak coupling 
to the leads the full retarded Green function is close to 
the unperturbed one and the electron dynamics inside 
the system must resemble the one of the isolated sample. 
Indeed, the curve at U = 0.25 follows the oscillations of 
the free Green function which is also given in the figure 
(the curve corresponding to U = 0.0). We plot the imag- 
inary part since it turns out again that the real part is 
vanishingly small so it does not contribute considerably 
to the retarded current. As U increases the Green func- 
tion changes and shows clear oscillations imposed by the 
leads' self-energy. 

The analysis performed so far was focused on the be- 
havior of the transient current as the intrinsic parameters 
of the system (i.e. its size and the height of the tunneling 
barriers at the contacts are varied. Nevertheless, in a typ- 
ical transport experiment these parameters are fixed and 
one usually measures the current by varying the bias or 
a plunger gate voltage. From steady-state current mea- 
surements it is well known that the role of such a gate 
potential is to bring one or more levels of the quantum 
dot within the bias window (BW) . We show in what fol- 
lows that at fixed bias and given coupling strength to the 
leads one can tune the transient current with a gate po- 



tential. Moreover, by inspecting the transient behavior 
as the gate potential is varied, it is possible to extract 
some information about the number of states within the 
bias window or above it. We will make the discussion 
for the 4 site dot. The gate potential is simulated by 
the diagonal term V g added to the on-site energy of the 
system. We fix the bias window to W=2.0 and for con- 
venience we set fin, = 0.0. Fig. 8(a) give a families of 
transients for coupling strength U = 0.75 and various 
values of V g specified in the figures. Fig. 8(b) shows the 
four levels of the isolated quantum dot as the gate po- 
tential scans the range [—4 : 4]. The V g values chosen 
in Fig. 8(a) corresponds to different location of the levels 
w.r.t. to BW. The bottom curve is irregular and settles 
down to a vanishing current because in this case there is 
no level within the BW. We note however that a nonva- 
nishing transient current still develops shortly after the 
coupling is established. At V g = —1.0 the highest level 
is located in the BW and the transient is smooth and al- 
ready shows the additional shoulder noticed previously 
The same thing happens when two states lie in the BW 
(at V g = 0.0), the difference being that the steady state 
current increases considerably. For V g > 0.50 it is clear 
from the structure of the spectrum that one cannot have 
more than two states in the BW and that the levels pass 
gradually above it. We found interesting to look at the 
transient currents for those gate potentials that still allow 
two states in the transmission range while pushing one 
or two states above BW. One notices that for V g = 1.0 
the steady state currents do not distinguish the different 
spectral structure involved in transport, while the tran- 
sient current is very sensitive to it. 

In the case of the six site QD the shoulder in Fig. 4 
is more pronounced because at V g = 0.0 there are ex- 
actly three states inside the bias window. We want to 
point out that since we have neglected the Coulomb in- 
teraction our simulations cannot capture the transport 
through excited states of the quantum dot. Tunneling 
processes involving such states would lead to a minipeak 
structure of the current maxima. We can consider how- 
ever our results should describe qualitatively the trans- 
port involving more levels because for small dots the bias 
required to cover the ground state of N electrons is much 
higher that the excitation energies. 

Now we investigate two more features of the transient 
regime: The time-dependent charge filling of the central 
region and possible effects due to different shapes of the 
system or of the various ways in which one can couple 
the leads. Besides the 4 site ID system discussed so far 
we consider also a 2 x 2 quantum dot. Both systems 
are submitted to the same bias and have equal coupling 
to the leads. However, in the case of the 2D quantum 
dot one can use different contacts for plugging the leads. 
We discuss two situations: i) A symmetric configuration 
in which the leads are attached to the opposite corners 
of the system, namely at the sites 1 and 4 and ii) An 
asymmetric coupling when we use the 1st and 3th sites 
as contacts. Fig. 9(a) reveals the changes induced in the 
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FIG. 8: (Color online) (a) The transient current for different 
values of the gate potential applied on the 4 site quantum dot 
for U — 0.75. (b) The spectrum of the system as function 
of the gate potential. It can be checked that the additional 
shoulders observed in (a) develop when there is at least one 
state of the QD in the bias window and no state above it. 
Other parameters: V = 2, and kT = 0.0001. Note that the 
bias is applied asymmetrically, that is up = 0.0. 
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transient curves in each case. The remaining subfigures 
show the occupation number Ni(t) of each site i = 1,4 
for the ID system and the two configurations considered. 
Inspecting Fig. 9(a) and (b) we see: i) The contact sites 1 
and 4 are the first to be populated due to their proximity 
to leads; ii) Since fi a > fip the right contact site (the 4- 
th) gains less charge at a lower rate than the left contact 
(the first). Both N\ and N4 show a step-like behavior 
for a short period (around t = 1). This coincides with 
the increased occupation number on the middle sites. We 
note also that the step in the occupation of the contact 
site N\ ends when it is equaled by N2- All the sites are 
then continuously filled up to the steady-state value. The 
occupation number on the right contact is smaller than 
the other ones which attain roughly the same value 0.65. 
iii) The step-like behavior of the transient currents in the 
range [1 : 1.75] corresponds to the almost constant popu- 
lation of the contact sites in the same interval. The sym- 
metric configurations is still characterized by a smooth 
transient but we notice that the step appears now later 




FIG. 9: (Color online) (a) The transient in the left lead J a (t) 
for the 4x1 system and for the 2x2 system in the two 
configurations of the leads as mentioned in the text, (b) - 
ID system, (c) - symmetric configuration, (d) - asymmetric 
The occupation numbers Ni(t) of the i-th site for the three 
cases considered in Fig. 8(a). Other parameters: U = 0.75, 
V = 2.0, and kT = 0.0001. 
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FIG. 10: (Color online) The total transient current and the 
components in the left lead and right lead. Other parameters: 
U = 0.5, V = 2, kT = 0.0001. 



that in the ID case. Fig. 9(c) confirms again that this sta- 
ble regime is assigned to a constant flow in the contact 
site. Also it reflects that the charge is equally distributed 
in the sites 2 and 3 which are located symmetrically w. 
r.t. the leads. In the asymmetric geometry the transient 
is rather similar to the one of the ID system up to t = 2.8 
but then drops to a lower steady state value. The occu- 
pation numbers show that the fourth site carries more 
charge than the contact sites in the steady state. This 
means in our opinion that part of this charge simply ac- 
cumulates and is not participating in transport. More 
interestingly, we note that in contrast to the symmetric 
geometry N% and N3 are different in the transient regime 
but reach the same value in the steady state. We mention 
that time-dependent simulations were performed recently 
in the case of an Aharonov-Bohm ring starting from the 
Schrodinger equationii^i 

We discuss now briefly the total current J a + Jp which 
is given in Fig. 10 along with its two components (the 
+ sign is due to the fact that the current Jp represents 
the current from the lead /3 to the system and therefore 
has opposite sign). As already mentioned, only in the 
steady state the current conservation reads as J a = — J p. 
Consequently, the net current J a + Jp vanishes. However, 
in the transient regime the two currents, although having 
similar shape differ significantly. 

To substantiate further the previous analysis we 
present in Fig. 11(a) the contour plot of the current as 
function of time and gate potential for the 2x2 site 
quantum dot coupled to the leads in the symmetrical 
configuration. Fig. 11(c) gives the spectrum of the sys- 
tem as the gate potential varies. The middle eigenvalue 
is doubly degenerated. When the levels are either be- 
low or above the bias window (W — 1.0 starting from 
jtiR = 0.0) the transient oscillates for quite a long time 
before passing to the steady state. We also observe that 
in these two extreme limits the transient oscillations are 
qualitatively different. For V g ~ —4 the current shows 
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FIG. 11: (Color online) (a) The 3D map of the transient 
current for a 2 x 2 site QD as a function of time and gate 
potential V g . The maxima are related to the spectrum of 
the isolated system given in (c). (b) The 3D map of the 
transient current for 2x2 the double dot t24 = £13 = 0.5 
Other parameters: V = 1.0, U = 0.50, kT = 0.0001. (c) The 
spectra of the two systems as a function of the gate potential. 
The lines mark the bias window. 



decreasing oscillations towards the steady state. In con- 
trast, for V g ~ 4 one remarks faster oscillations and more 
importantly, negative values of the transient. Since in 
this regime all the levels are above the bias window and 
there is no way to pass electrons to the right lead it is 
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clear that the negative current in the left lead is just the 
reflected one. We underline that this effect is due to the 
fact that we have considered a finite spectral width of 
the leads and that similar features were reported for the 
single site case^. As expected, as the system approaches 
the stationary regime the current shows three maxima 
associated to the passage of the localized levels through 
the bias window. Actually the levels turn to resonances 
when coupling the leads to the system but since one has 
to deal with a time-dependent Hamiltonian it is difficult 
to characterize the location and width of these resonances 
in the transient regime. This is why in the 3D plot one 
cannot distinguish between different resonances at times 
t < 3.0. Fig. 11(b) presents the transient current for the 
same 2x2 system except that the hopping parameters £13 
and ^24 are reduced to 0.5. In this case one can view the 
system as a double dot, each dot composed of two sites. 
As the spectrum from Fig. 11(b) shows, the degeneracy is 
lifted and the level spacing diminishes. As a consequence 
in the long time regime one gets two broader peaks, since 
the four levels are now grouped into pairs. 

All the features presented above emphasize that the 
transient regime of the many-level structures is quite dif- 
ferent from the single-level system. 

IV. CONCLUSIONS 

We have performed transient current calculations for 
a many-level finite system coupled suddenly to semiinfi- 
nite biased leads. Our method is based on the nonequi- 
librium Green-Keldysh machinery. We find numerically 



an exact solution of the integral Dyson equation which is 
solved as an algebraic equation. By analyzing the behav- 
ior of the retarded and lesser Green functions we explain 
qualitatively the shape of the transient current and the 
passage to the steady state. The amplitude of the cou- 
pling to the leads controls essentially the convergence to 
a steady state. We have identified non-trivial effects of 
the many-level structure of the system and presented an 
intuitive picture of the charge filling by studying the oc- 
cupation number inside the system. By increasing the 
system size the shape of the transient current and the 
evolution towards the steady state differs significantly 
from the single-site oscillatory behavior and depends cru- 
cially on the number of electronic states available in the 
bias window. We predict that a step-like structure could 
be observed in transient current measurements by apply- 
ing a gate potential on the system that tunes the higher 
levels within the bias window. Different transients are 
expected to appear as well when different coupling ge- 
ometries of the leads are used. 

The present method can be used for studying the re- 
sponse of mesoscopic systems to more complicated time- 
dependent couplings to the leads: pulses having different 
lengths and decaying rates, and non-periodic signals. 
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